#clean up
rm(list=ls())

#Set Working Directory
#setwd('/Volumes/MONOGAN/stateAreal/')
#setwd("C:/Users/Josh Jackson/Dropbox/Spatial Paper with Jamie/Replication")

####Decision Tree Graph####
#postscript("decisionTree.eps",horizontal=FALSE, width=5, height=5, paper='special', onefile=FALSE, pointsize=10)
#pdf("decisionTree.pdf",width=5, height=5, pointsize=12)
plot(x=c(0,1.25),y=c(0,1),xlab="",ylab="",axes=F,type="n")
arrows(x0=.5,y0=.9,x1=.25,y1=.5,length=.1)
arrows(x0=.5,y0=.9,x1=.75,y1=.5,length=.1)
text(x=.5,y=.9,label="Functional Form",pos=3)
text(x=.35,y=.75,label="Diffusive",pos=2)
text(x=.65,y=.75,label="Confined",pos=4)
text(x=.25,y=.5,label="Spatial Lag Model",pos=1)
text(x=.775,y=.5,label="Error Structure",pos=1)
arrows(x0=.75,y0=.4,x1=.6,y1=.1,length=.1)
arrows(x0=.75,y0=.4,x1=.9,y1=.1,length=.1)
text(x=.675,y=.25,label="Direct",pos=2)
text(x=.825,y=.25,label="Hierarchical",pos=4)
text(x=.5,y=.1,label="MLE CAR/SAR",pos=1)
text(x=1,y=.1,label="Bayesian CAR",pos=1)
#dev.off()


####Forest Plot Results####
#Cross-Sectional Results
#EWM RESULTS:
#postscript("ewmForest.eps",horizontal=FALSE, width=5, height=5, paper='special', onefile=FALSE, pointsize=10)
#pdf("ewmForest.pdf",width=5, height=5,  pointsize=10)
lin.mean<-c(.7626, .8414, -1.023, .1294, -.2245, .7866, 1.053, .6601, .4428, -.3308, .4775)
lin.lower<-c(.5964, .5443, -1.32, -.1139, -.4867, .6171, .97, .5772, .2218, -.4606, .2508)
lin.upper<-c(.9294, 1.138, -.7259, .372, .03814, .9564, 1.136, .7431, .6647, -.2013, .7049)

car.mean<-c(0.6033, 0.8221,-0.8485,0.06182,-0.1247,0.6976,1.066,0.6652,0.3687,-0.2959,0.4933)
car.lower<-c(0.4052,0.4937,-1.215,-0.2279,-0.4378,0.5007,0.9106,0.5215,0.09225,-0.4948,0.2178)
car.upper<-c(0.8015,1.152,-0.4791,0.3506,0.1898,0.8948,1.221,0.8088,0.6458,-0.0954,0.7687)

abs(car.upper-car.lower)-abs(lin.upper-lin.lower)
   #CAR credible intervals wider across the board

par(omi=c(.2,1.1,.1,.01))
plot(x=lin.mean,y=c(1:11-.1),xlim=c(-1.3,1.3),ylim=c(0,12), xlab="Coefficient", ylab="",axes=F)
axis(1)
axis(2, at=c(1:11), labels=c('Opinion on party elite','Opinion on Dem ID','Party elite on Dem ID','Opinion on strength','Party elite on strength','Dem ID on strength','Party elite on legis lib','Strength on legis lib','Opinion on policy','Strength on policy','Legis lib on policy'),las=1)
box()
abline(v=0,col='gray60')
segments(x0=lin.lower,x1=lin.upper,y0=c(1:11)-.1,y1=c(1:11)-.1, lty=2)
points(x=car.mean,y=c(1:11+.1), col='red',pch=2)
segments(x0=car.lower,x1=car.upper,y0=c(1:11+.1),y1=c(1:11+.1), col='red')
legend(x=-1.3, y=2.5, pch=c(2,1), lty=c(1,2), col=c('red', 'black'), legend=c('CAR', 'Independent'))
#dev.off()

#Duration Results
#Lottery RESULTS:
#postscript("lotteryForest.eps",horizontal=FALSE, width=5, height=5, paper='special', onefile=FALSE, pointsize=10)
#pdf("lotteryForest.pdf",width=5, height=5,  pointsize=10)
lin.mean<-c(.6096,.504,-.004703,.1593,-.1271,.5193,.05508,-.02823)
lin.lower<-c(-.1845,-.224,-.1197,-.4473,-.1786,.2399,.008679,-.04531)
lin.upper<-c(1.43,1.279,.1307,.7681,-.0763,.7999,.09814,-.0120)

car.mean<-c(0.5633,.5017,.05112,0.1383,-.103,0.08812,.06677,-.03728)
car.lower<-c(-.2246,-.2354,-.0601,-.4904,-.1584,-.2188,.0249,-.05607)
car.upper<-c(1.388,1.28,.1593,0.7553,-.05092,0.3978,.1095,-.01963)

par(omi=c(.2,.9,.1,.01))
plot(x=lin.mean,y=c(1:8-.1),xlim=c(-.75,1.5),ylim=c(0,11), xlab="Coefficient", ylab="",axes=F)
axis(1)
axis(2, at=c(1:8,10), labels=c('Election 1','Election 2','Lagged income','Unified party','Fundamentalists','Neighbors adopted','Time','Time squared', 'Lagged fiscal health'),las=1)
box()
#abline(v=0,col='gray60')
segments(x0=0,x1=0,y0=-1,y1=8.6, col='gray60')
segments(x0=lin.lower,x1=lin.upper,y0=c(1:8-.1),y1=c(1:8-.1), lty=2)
points(x=car.mean,y=c(1:8+.1), col='red',pch=2)
segments(x0=car.lower,x1=car.upper,y0=c(1:8+.1),y1=c(1:8+.1), col='red')
#legend(x=-1, y=1, pch=c(2,1), lty=c(2,1), col=c('red','black'), legend=c('CAR', 'Independent'))
abline(h=8.6,col='black')
abline(h=8.7,col='black')

abs(car.upper-car.lower)-abs(lin.upper-lin.lower)
#mean(abs(car.upper-car.lower)-abs(lin.upper-lin.lower))
   #CAR wider on 5 out of 8. 

par(new=T,omi=c(.2,.9,.1,.01))
lin.mean<-c(-3.422)
lin.lower<-c(-6.847)
lin.upper<-c(-.2447)

car.mean<-c(-3.343)
car.lower<-c(-6.817)
car.upper<-c(-.08333)

abs(car.upper-car.lower)-abs(lin.upper-lin.lower)
   #inconsistent
	#CAR wider on 6 out of 9.

plot(x=lin.mean,y=c(10-.1),xlim=c(-9,4),ylim=c(0,11), xlab="Coefficient", ylab="",axes=F)
axis(3)
segments(x0=0,x1=0,y0=8.7,y1=12, col='gray60')
segments(x0=lin.lower,x1=lin.upper,y0=c(10-.1),y1=c(10-.1), lty=2)
points(x=car.mean,y=c(10+.1), col='red',pch=2)
segments(x0=car.lower,x1=car.upper,y0=c(10+.1),y1=c(10+.1), col='red')
legend(x=-1.75, y=8.5, pch=c(2,1), lty=c(1,2), col=c('red','black'), legend=c('CAR', 'Independent'))
#dev.off()


#Multilevel model results (Margalit 2013) 
#postscript("margalitForest.eps", width=5, height=5, paper='special', onefile=F, pointsize=10)
#pdf("margalitForest.pdf", width=5, height=5, pointsize=10)
sw.mean=   c(.041,-.006, .000, .018, .171,-.122, .003,-.005, .009,-.008)
sw.lower=  c(.024,-.022,-.017, .003, .146,-.145,-.013,-.020,-.021,-.025)  
sw.upper=  c(.060, .010, .017, .033, .195,-.098, .019, .009, .038, .008)

car.mean=  c(.042,-.006, .000, .018, .171,-.122, .003,-.005, .009,-.006)
car.lower= c(.024,-.022,-.018, .003, .146,-.145,-.013,-.020,-.021,-.021)
car.upper= c(.060, .010, .018, .033, .195,-.098, .019, .010, .038, .010)

par(omi=c(.1,.8,.1,.1))
plot(x=sw.mean, y=c(1:10)-.15, ylim=c(0,11.5), xlim=c(-.15,.2), xlab="Coefficient", ylab="", axes=F)
axis(1,cex.axis=.8)
axis(2, at=c(1:11), labels=c("Lost Job", "HH Income Drop", "Job Less Secure", "Sp. Lost Job", "Democrat", "Republican", "Long Term Unemp.", "New Job", "Not Looking", "Unemployment", "Prev. Welfare"), las=1)
box()
#abline(v=0, col="gray60")
segments(x0=sw.lower,x1=sw.upper,y0=c(1:10)-.15, y1=c(1:10)-.15, lty=2)
points(x=car.mean,y=c(1:10)+.15, col='red',pch=2)
segments(x0=car.lower,x1=car.upper,y0=c(1:10)+.15,y1=c(1:10)+.15, col='red')
abline(h=10.6, col="black")
abline(h=10.5, col="black")
segments(x0=0,x1=0,y0=-1,y1=10.5, col='gray60')
legend(x=0.05, y=9, pch=c(2,1), lty=c(1,2), col=c('red', 'black'), legend=c('CAR', 'Independent'))

#"Prev. Welfare"
par(new=T, omi=c(.1,.8,.1,.1))
sw.mean2= .588
sw.lower2= .569
sw.upper2= .608

car.mean2= .589
car.lower2= .569
car.upper2= .608

plot(x=sw.mean2,y=c(11-.15),ylim=c(0,11.5),xlim=c(0,0.65), xlab="Coefficient", ylab="",axes=F)
axis(3, cex.axis=.8)
segments(x0=0,x1=0,y0=10.6,y1=12, col='gray60')
segments(x0=sw.lower2,x1=sw.upper2,y0=c(11-.15),y1=c(11-.15), lty=2)
points(x=car.mean2,y=c(11+.15), col='red',pch=2)
segments(x0=car.lower2,x1=car.upper2,y0=c(11+.15),y1=c(11+.15), col='red')
#dev.off()

abs(car.upper-car.lower)-abs(sw.upper-sw.lower)
   #The CAR interval is wider in two cases, and the independent is wider in one. In all other cases, the intervals are equivalent to our level of precision.
abs(car.upper2-car.lower2)-abs(sw.upper2-sw.lower2)
   #Identical 

#Clean Air results  
#postscript("cleanairForest.eps", width=5, height=5, paper='special', onefile=F, pointsize=10, family='ComputerModern')
#pdf("cleanairForest.pdf", width=5, height=5, pointsize=10)

ind.mean= c(.0709, .0572, .0222, .0040, .1836, .0278, .0014)
ind.lower=c(.0624, .0502, .0145,-.0023, .1667, .0230,-.0046)
ind.upper=c(.0792, .0640, .0299, .0101, .1999, .0325, .0073)

car.mean= c(.0705, .0573, .0219, .0038, .1858, .0278, .0015)
car.lower=c(.0604, .0488, .0127,-.0037, .1642, .0219,-.0058)
car.upper=c(.0804, .0657, .0313, .0112, .2069, .0335, .0089)

par(omi=c(.1,.8,.1,.1))
plot(x=ind.mean, y=c(1:7)-.15, ylim=c(0,7.5), xlim=c(-.015, .215), xlab="Coefficient", ylab="", axes=F)
axis(1)
axis(2, at=c(1:7), labels=c("Fed. Actions (lag)", "Dems Unified", "Reps Unified", "Unemployment", "Manufacturing", "Grants", "Groups"), las=1)
box()
segments(x0=ind.lower,x1=ind.upper,y0=c(1:7)-.15, y1=c(1:7)-.15, lty=2)
points(x=car.mean,y=c(1:7)+.15, col='red',pch=2)
segments(x0=car.lower,x1=car.upper,y0=c(1:7)+.15,y1=c(1:7)+.15, col='red')
abline(v=0, col='gray60')
legend(x=0.12, y=1, pch=c(2,1), lty=c(1,2), col=c('red', 'black'), legend=c('CAR', 'Independent'))
#dev.off()

abs(car.upper-car.lower)-abs(ind.upper-ind.lower)
#wider across the board

####RAW DATA PLOTS####
#Lottery Data
#setwd("/Volumes/MONOGAN/stateAreal/lottery/")
bbNew<-read.csv("bbExtended4.csv")
bbNew$zInc<-(bbNew$lagIncome-mean(bbNew$lagIncome))/sd(bbNew$lagIncome)
bbNew$zFisc<-(bbNew$lagFiscal-mean(bbNew$lagFiscal))/sd(bbNew$lagFiscal)
bbNew$zRelig<-(bbNew$religion-mean(bbNew$religion))/sd(bbNew$religion)
lotto<-bbNew
table(lotto$year)
count<-abs(table(lotto$year)-50)[-1]
count<-c(0,count,43,43,43,43,44)
years<-c(1963:2017)
names(count)<-years

#postscript("lotto.eps",horizontal=FALSE, width=5, height=5, paper='special', onefile=FALSE, pointsize=10)
#pdf("lotto.pdf",width=5, height=5, pointsize=10)
plot(y=count,x=years,xlab="Year",ylab="Number of States With Lotteries",type='o',ylim=c(0,50),lwd=2,pch=20)
abline(h=0,col='gray60')
#dev.off()


#Clean Air Enforcement Data
library(lattice)
library(foreign)
#setwd("/Volumes/MONOGAN/stateAreal/wood1992/")
caa<-read.dta('cleanAir0109.dta')
caa<-subset(caa, year>=2001)
caa<-subset(caa, year<=2009)

#trellis.device("pdf",color=FALSE,file="stateActions0109.pdf")
#trellis.device("pdf",color=TRUE,file="stateActions0109.pdf")
#postscript("stateActions.eps", horizontal=FALSE, colormodel='rgb', width=5, height=5, paper="special", onefile=FALSE, family='ComputerModern')
xyplot(log(stateact)~year, caa, type='l', groups=state, xlab="Year", ylab="State Actions (Logged Scale)")
#dev.off()
